% function BuildTimingDiagram()
%BUILDTIMINGDIAGRAM
logfiles = {};
logfiles{end+1} = 'D:\Ralph\ECoG Code\Playground\BuildTimingDiagram.txt';
diary('D:\Ralph\ECoG Code\Playground\BuildTimingDiagram.txt');

% Approach:
%  V Get the preprocessed data for subjects
%  V For all conditions and regions, merge the data over all subjects (Create a new EEGLAB set for this data!)
%  V Epoch the data

%   (This is a single trial per condition per region)
%   TFC this data and do peakfinding
%   Plot the diagram

Datapath = [getpref('ECoG', 'Data') 'ByRegion\'];
output.datatable = [];

Subjects;
% TO LIMIT THE SUBJECTLIST, ADD VARIABLE subjectselector = [..] before
% running the script
fn = fieldnames(subject);
if ~exist('subjectselector', 'var'), subjectselector = 1:length(fn); end;
for i=1:length(fn(subjectselector))
    subjectinfo = subject.(fn{subjectselector(i)});
    
    %% Generate merged region data for each subject
    datasetFilenames = {};
    for d=1:length(subjectinfo.Datasets) % for every dataset
        EEG = pop_loadset('filename', subjectinfo.Datasets(d).PreprocessedFilename);
        EEG = eeg_checkset( EEG );
        [REGEEG ChannelInfo] = MergeEEGDATAByRegion(subjectinfo, EEG);
        diary off;
        logfiles{end+1} = 'D:\Ralph\ECoG Code\Playground\Result of MergeEEGData.txt';
        diary('D:\Ralph\ECoG Code\Playground\Result of MergeEEGData.txt');
        disp('--------------------------------------------------------------------');
        disp(['Resulting merges for ' subjectinfo.Name ' in dataset ' subjectinfo.Datasets(d).PreprocessedFilename ':']);
        for chIdx = 1:length(ChannelInfo)
            disp([' -- ' ChannelInfo(chIdx).Name ', ' num2str(ChannelInfo(chIdx).ElectrodeCount) ' electrodes']);
            eleclist = '';
            for elecI = 1:length(ChannelInfo(chIdx).Electrodes), eleclist = [eleclist ' ' ChannelInfo(chIdx).Electrodes(elecI).Label]; end;
            disp(['        ' eleclist]);
        end
        diary off;
        diary('D:\Ralph\ECoG Code\Playground\BuildTimingDiagram.txt');
        datasetFilenames{d} = [Datapath 'Reg' subjectinfo.PrimaryPrefix ' Set ' num2str(d) '.set'];
        pop_saveset(REGEEG, 'filename', datasetFilenames{d});
        EEG = []; REGEEG = []; % Todo, improve this, it is fragmenting the memoryspace..
    end % end for every dataset
    cfgBE = buildEpochs_defaultCfg();
    cfgBE.DatasetFilenames = datasetFilenames'; % Instead of the EEGlab datasets in subjectinfo, use these filenames
    cfgBE.Output.AutogenerateFilenames = false; % Set this to false to use the prefix and directory below
    cfgBE.Output.EpochsetFilenamePrefix = ['Reg' subjectinfo.PrimaryPrefix]; % Prefix the output filename with this.
    cfgBE.Output.EpochsetDirectory = Datapath; % Prefix the output filename with this.
    cfgBE.Output.StoreInSubjectinfo = false; % When true, buildEpochs will rewrite the EpochData node in the subjectinfo
    cfgBE.SaveSubjectinfo = false;
    [outputBE subjectinfo] = buildEpochs(subjectinfo, cfgBE);
    diary off;
    logfiles{end+1} = 'D:\Ralph\ECoG Code\Playground\Result of BuildEpochs.txt';
    diary('D:\Ralph\ECoG Code\Playground\Result of BuildEpochs.txt');
    disp('--------------------------------------------------------------------');
    disp(['Resulting epochs for subject: ' subjectinfo.Name]);
    fnEP = fieldnames(outputBE.Datasets);
    for fnEPi = 1:length(fnEP)
        disp([' -- ' num2str(outputBE.Datasets.(fnEP{fnEPi}).nrOfEpochsFound) ' epochs for ' outputBE.Datasets.(fnEP{fnEPi}).name]);
    end
    diary off;
    diary('D:\Ralph\ECoG Code\Playground\BuildTimingDiagram.txt');
    
    datasetsTotal = fieldnames(outputBE.Datasets)';
    disp(datasetsTotal');
%     if ~exist('datasetsToDoSelector', 'var'),
        datasetsToDoSelector = 1:length(datasetsTotal); 
%     end;
    disp(datasetsToDoSelector);
    % USE ABOVE TO TEMPORARILY SELECT A SUBSET OF DATA TO DO..
    for fnBE = datasetsTotal(datasetsToDoSelector)   % run over conditions (=epochsets)
        x = outputBE.Datasets.(cell2mat(fnBE));
        disp(x.name);
        
        %% Generate timefreq data for each epochset
        cfg = buildFrequencySet_defaultCfg();
        cfg.SubjectPrimaryPrefix = ['Reg' subjectinfo.PrimaryPrefix];
        cfg.EpochSetName = x.name;
        cfg.ProcessAllChannels = false;
        cfg.ChannelSelection = {ChannelInfo.Name}; % Use this to select which channels to process in case ProcessAllChannels is false
        cfg.ChannelSelection = cfg.ChannelSelection ([ChannelInfo.ElectrodeCount] > 0);
        cfg.Input.EegSetFilename = x.dataFilename;
        cfg.Output.TfcDataFolder = [getpref('ECoG','Data') 'ByRegion\Timefreq\'];
        disp(cfg);
        outputBFS = buildFrequencySet_Core(cfg);
        
        %% Generate tfc data for each epochset
        fnCh = fieldnames(outputBFS.TfcData);
        listOfTfcDataFiles = cell(1,length(fnCh));
        for j = 1:length(fnCh) % for every channel
            listOfTfcDataFiles{j} = outputBFS.TfcData.(fnCh{j}).Filename;
        end % end for every channel
        cfg = calculateTfcData_defaultCfg();
        cfg.SubjectPrimaryPrefix = ['Reg' subjectinfo.PrimaryPrefix];
        cfg.SubjectName = subjectinfo.Name;
        cfg.EpochSetName = x.name;
        cfg.Input.ListOfTfcFiles = listOfTfcDataFiles;
        cfg.Output.TfcDataFolder = [getpref('ECoG','Data') 'ByRegion\Timefreq\'];
        cfg.RemoveEpochsConditionally = true;
        disp(cfg);
        diary off;
        logfiles{end+1} = 'D:\Ralph\ECoG Code\Playground\Results of CalculateTfcData.txt';
        diary('D:\Ralph\ECoG Code\Playground\Results of CalculateTfcData.txt');
        outputCTD = calculateTfcData_Core(cfg);
        diary off;
        diary('D:\Ralph\ECoG Code\Playground\BuildTimingDiagram.txt');
        
        %% Generate peak information
        currentEpochset = x.taskname;
        nrOfEpochs = x.nrOfEpochsFound;
        TfcDatafile = outputCTD.TfcDatafile;
        % Load the data
        load(TfcDatafile);
        % Set the configuration
        % cfg = scoreTfcPeaks_defaultCfg();
        % Score the peaks
        channels = fieldnames(tfc.data);
        warning('off', 'all');
        for gpiCh = 1:length(channels) % for every electrode
            channel = channels{gpiCh};
            for gpiFrq=1:length(tfc.data.(channel)) % for every freq.bin
                PEAKS = scorePeaks(tfc.data.(channel)(gpiFrq), tfc.times, nrOfEpochs);
                tfc.data.(channel)(gpiFrq).peaklocs = PEAKS.peaklocs;
                tfc.data.(channel)(gpiFrq).peaktype = PEAKS.peaktype;
                tfc.data.(channel)(gpiFrq).firstEffect = PEAKS.firstEffect;
                tfc.data.(channel)(gpiFrq).pvalue = PEAKS.pvalue;
            end % end for every freq.bin
        end % end for every electrode
        warning('on', 'all');
        % Update the datafile
        save(TfcDatafile, '-v7.3', 'tfc');
        
        %% Rewrite the tfcdata peaks into a datatable
        [datatable header] = TfcToDatatable(tfc);
        
        if isempty(datatable), output.datatable = datatable;
        else output.datatable = [output.datatable; datatable];
        end
        
        % Save stuff for now
        save(['Temporary ' subjectinfo.PrimaryPrefix x.name '.mat'], '-v7.3')
        
    end % end for every condition
    
end % end for every subject


%% Plot all timings
% bins = {'4-7Hz', '7-13Hz', '13-30Hz', '30-70Hz', '70-130Hz'};
% pptCfg = plotPeakTiming_defaultCfg();
% for b = 1:length(bins)
%     figure;
%     pptCfg.binSelection = bins{b};
%     for i = 1:length(setToDo)
%         epochset = setToDo{i};
%         load([epochset '.mat'], 'datatable', 'header');
%         subplot(length(setToDo),1,i);
%         plotPeakTiming(datatable, header, pptCfg);
%         title(epochset);
%     end
%     a = axes('position',[.02 .02 .96 .96],'Visible', 'off');
%     text('Position', [0.01 0.99], 'String', [subjectinfo.Name ' / ' bins{b}], ...
%         'HorizontalAlignment', 'left', 'VerticalAlignment', 'cap', 'Parent', a);
%     text('Position', [0.5 0.01], 'String', datestr(now, 'yyyy-mm-dd HH:MM:SS'), ...
%         'HorizontalAlignment', 'center', 'VerticalAlignment', 'baseline', 'Parent', a);
% end;

% end

setpref('Internet','SMTP_Server','smtp.uu.nl');
setpref('Internet','E_mail','r.j.m.hendriks@uu.nl');
sendmail({'ralph@affy.nl','r.j.m.hendriks@uu.nl','dejongmc@gmail.com'},'Matlab diary','See attachement for transscript', logfiles);
